#define AnaTask1_cxx
#include "AnaTask1.h"

#include "ND280HL2Utils.h"

void AnaTask1::DoAnalysis(const TString filelist, const Int_t ntarget, const TString tag, const Int_t krec)
{
  NeutrinoUtils::IniCCQE(); //Initilises all the parameters for NeutrinoUtils - called from utils

  Init(ND280HL2Utils::IniIO(filelist, krec?"default":"truth")); //Sets all the branch pointers and initialises called from header AnaTask1.h

  printf("AnaTask1::DoAnalysis ntarget %d\n", ntarget);

  TList *lout=NeutrinoUtils::GetHistList(); // Sets the list of callable histgrams

  TTree *tout = NeutrinoUtils::GetTree(); // Sets branch location to TTrees

  const Int_t nentries = fChain->GetEntries();//fChain->GetEntriesFast(); - where is this called from check .h file
  Int_t nsel=0;
  Int_t nfail = 0;
  for (Int_t jentry=0; jentry<nentries;jentry++) {
    if(jentry%10000==0){
      printf("%d/%d nsel %d nfail %d\n", jentry, nentries,nsel, nfail);
    }

    if(!fChain->GetEntry(jentry)){
      printf("AnaTask1::DoAnalysis GetEntry ends at jentry %d\n", jentry);
      break;
    }

    //should be similar to numuCCAnalysis
    if(accum_level[0][1]<11){
      nfail++;
      continue;
    }

    //--------------------------------------------
    NeutrinoUtils::FillCount(lout, 0);

    ND280HL2Utils::SaveInfo(nu_truereac, nu_parentpdg, nu_pdg, ntpctracks);

    //---------------------------------------------------------------------------------------------------------
    //reconstruction
    //---------------------------------------------------------------------------------------------------------

    ND280HL2Utils::SetNeutrinoRec(vertex_pos); // Sets the Neutrino interaction point returning NeutrinoUtils::fNeutrinoRec

//ND280HL2Utils::SetParticleRec(NeutrinoUtils::fNeutrinoRec, NeutrinoUtils::fMuonRec,   NeutrinoTools::MuonMass(),   particle direction [3], mom magnitude, vertex_pos, selmu_endpos, selmu_Chi2[0], selmu_NDOF[0], selmu_ntpcs, selmu_nfgds, selmu_tpc_nhits, selmu_tpc_charge);

    const Int_t nmurec = ND280HL2Utils::SetParticleRec(NeutrinoUtils::fNeutrinoRec, NeutrinoUtils::fMuonRec,   NeutrinoTools::MuonMass(),   selmu_dir, selmu_mom[0], vertex_pos, selmu_endpos, selmu_Chi2[0], selmu_NDOF[0], selmu_ntpcs, selmu_nfgds, selmu_tpc_nhits, selmu_tpc_charge);
    NeutrinoUtils::SetMuonNRec(nmurec);

    const Int_t nprec  = ND280HL2Utils::SetParticleRec(NeutrinoUtils::fNeutrinoRec, NeutrinoUtils::fProtonRec, NeutrinoTools::ProtonMass(), selp_dir,  selp_mom[0],  vertex_pos, selp_endpos,  selp_Chi2[0],  selp_NDOF[0],  selp_ntpcs,  selp_nfgds,  selp_tpc_nhits,  selp_tpc_charge);

    NeutrinoUtils::SetProtonNRec(nprec);

    //---------------------------------------------------------------------------------------------------------
    //simulation 
    //---------------------------------------------------------------------------------------------------------

    ND280HL2Utils::SetNeutrinoSim(nu_truedir, nu_trueE);

    const Int_t nmusim = ND280HL2Utils::SetParticleSim(NeutrinoUtils::fNeutrinoSim, NeutrinoUtils::fMuonSim,   NeutrinoTools::MuonMass(),   truemu_truedir, truemu_truemom, particle);
      
    NeutrinoUtils::SetMuonNSim(nmusim);

    //proton: truep is the one out of nucleus
    const Int_t npsim  = ND280HL2Utils::SetParticleSim(NeutrinoUtils::fNeutrinoSim, NeutrinoUtils::fProtonSim, NeutrinoTools::ProtonMass(), truep_truedir,   truep_truemom,   selp_pdg);

    NeutrinoUtils::SetProtonNSim(npsim);

    //---------------------------------------------------------------------------------------------------------
    //Det
    //---------------------------------------------------------------------------------------------------------

    //proton: selp is the one after propagation in detector, suffering from ordinary multiple scattering
    ND280HL2Utils::SetParticleSim(NeutrinoUtils::fNeutrinoSim, NeutrinoUtils::fProtonDet, NeutrinoTools::ProtonMass(), selp_truedir,   selp_truemom,   selp_pdg);

    //---------------------------------------------------------------------------------------------------------
    //Toy
    //---------------------------------------------------------------------------------------------------------
	// Fermi Momtion toy
    ND280HL2Utils::SetProtonToy(0, NeutrinoUtils::fNeutronToy, NeutrinoUtils::fNeutrinoSim, NeutrinoUtils::fMuonSim->fLorentzVec, NeutrinoUtils::fProtonToy, selp_pdg);


    //---------------------------------------------------------------------------------------------------------
    //rec sim toy Delta
    //---------------------------------------------------------------------------------------------------------

    NeutrinoUtils::SetDeltaPt();
    NeutrinoUtils::SetpTT();

    //---------------------------------------------------------------------------------------------------------
    //---------------------------------------------------------------------------------------------------------

    NeutrinoUtils::FillCount(lout,11);

    tout->Fill();
    
    nsel++;

    if(ntarget>=0 && nsel>=ntarget){
      break;
    }
  }

  printf("nentries %d nsel %d nfail %d delta %d\n", nentries, nsel, nfail, nentries-nsel-nfail);

  if( (nentries-nsel-nfail)!=0 ){
    printf("AnaTask1 nsel nfail bad nsel %d nfail %d nentries %d\n", nsel, nfail, nentries);
    exit(1);
  }

  TFile *fout=new TFile(Form("ND280HL2%s_%d_%d_krec%d.root", tag.Data(), nsel, nfail, krec),"recreate");

  tout->Write();
  //for further handling in plot
  lout->Write("lout",TObject::kSingleKey);

  fout->Save();
  fout->Close();
  delete fout;
}

